drawing
Matthew Abbott / The New York Times


Introdução

Os incêndios florestais na Austrália são uma ocorrência regular que teve um papel significativo na formação da natureza do continente ao longo de milhões de anos. No entanto, os incêndios podem causar danos materiais significativos e levar à perda de vidas humanas e animais. Os incêndios florestais mataram cerca de 800 pessoas e bilhões de animais desde 1851. Estima-se que entre 2019 e 2020 os incêndios mataram pelo menos 33 pessoas e mais de 3 bilhões de animais.

Geralmente, os incêndios mais devastadores são precedidos por altas temperaturas, baixa umidade relativa e fortes ventos, que criam as condições ideais para a rápida propagação do fogo. A chuva é a maior aliada ao combate ao fogo e na mitigação do surgimento de novos de focos de incêndios. Outro problema que agrava a situação na Austrália, é a forte dependência do país de civis para conter seus incêndios, cerca de 90% dos bombeiros combatendo fogos florestais são voluntários.

Portanto, considerando os recursos limitados para combate aos incêndios florestais, ainda mais com o aumento das mudanças climáticas que estão tornando esses eventos mais frequentes e devastadores, saber se irá chover em uma região específica dadas algumas medições é importante para direcionar os esforços onde eles são mais necessários.

Visualizando os Dados

dataset = tibble(read.csv('weatherAUS.csv'))
print(dataset)
## # A tibble: 145,460 x 23
##    Date       Location MinTemp MaxTemp Rainfall Evaporation Sunshine WindGustDir
##    <chr>      <chr>      <dbl>   <dbl>    <dbl>       <dbl>    <dbl> <chr>      
##  1 2008-12-01 Albury      13.4    22.9      0.6          NA       NA W          
##  2 2008-12-02 Albury       7.4    25.1      0            NA       NA WNW        
##  3 2008-12-03 Albury      12.9    25.7      0            NA       NA WSW        
##  4 2008-12-04 Albury       9.2    28        0            NA       NA NE         
##  5 2008-12-05 Albury      17.5    32.3      1            NA       NA W          
##  6 2008-12-06 Albury      14.6    29.7      0.2          NA       NA WNW        
##  7 2008-12-07 Albury      14.3    25        0            NA       NA W          
##  8 2008-12-08 Albury       7.7    26.7      0            NA       NA W          
##  9 2008-12-09 Albury       9.7    31.9      0            NA       NA NNW        
## 10 2008-12-10 Albury      13.1    30.1      1.4          NA       NA W          
## # ... with 145,450 more rows, and 15 more variables: WindGustSpeed <int>,
## #   WindDir9am <chr>, WindDir3pm <chr>, WindSpeed9am <int>, WindSpeed3pm <int>,
## #   Humidity9am <int>, Humidity3pm <int>, Pressure9am <dbl>, Pressure3pm <dbl>,
## #   Cloud9am <int>, Cloud3pm <int>, Temp9am <dbl>, Temp3pm <dbl>,
## #   RainToday <chr>, RainTomorrow <chr>

Usando a API do Google Maps para salvar as coordenadas geográficas de cada local de medição, como existem vários lugares com o mesmo nome em diferentes países anglófonos, é necessário concatenar o nome do país no final do local para garantir que as coordenadas estarão coerentes

LocCount <- dataset %>% count(Location)

LocCount$Category <- paste(LocCount$Location, ", Australia")# Adding Australia to the end of location name since there are many cities with the same names in other countries

cities_df <- as.data.frame(LocCount)
locations_df <- mutate_geocode(cities_df, Category)
locations <- as_tibble(locations_df)

print(locations)
## # A tibble: 49 x 5
##    Location          n Category                    lon   lat
##    <chr>         <int> <chr>                     <dbl> <dbl>
##  1 Adelaide       3193 Adelaide , Australia       139. -34.9
##  2 Albany         3040 Albany , Australia         118. -35.0
##  3 Albury         3040 Albury , Australia         147. -36.1
##  4 AliceSprings   3040 AliceSprings , Australia   134. -23.7
##  5 BadgerysCreek  3009 BadgerysCreek , Australia  151. -33.9
##  6 Ballarat       3040 Ballarat , Australia       144. -37.6
##  7 Bendigo        3040 Bendigo , Australia        144. -36.8
##  8 Brisbane       3193 Brisbane , Australia       153. -27.5
##  9 Cairns         3040 Cairns , Australia         146. -16.9
## 10 Canberra       3436 Canberra , Australia       149. -35.3
## # ... with 39 more rows
df <- merge(x = dataset, y = locations, by = "Location", all = TRUE) # Outer join
# Creating Map Canvas
aus_coord <- geocode("Australia")

map <- get_googlemap(center = c(aus_coord$lon, aus_coord$lat), zoom = 4)

Mapa com as Temperaturas Máximas de Cada Local

MTemp <- df[!is.na(df$MaxTemp), ] %>% # Removing missing data
  select(Location, MaxTemp, lon, lat)


aggTemp <- aggregate(MTemp$MaxTemp, by = list(Locaton=MTemp$Location , Longitude=MTemp$lon , Latitude=MTemp$lat), FUN=mean) # Aggregating by Location and the Average  MaxTemp


ggmap(map) +
  geom_point(data = aggTemp,
             aes(x = Longitude, y = Latitude, size = x),
             color = "red", alpha = 0.3)  + 
  theme(legend.position = "left") +
  xlab("Longitude") + 
  ylab("Latitude") +
  labs(size="Average Max. Temperature (°C)")

Mapa de Calor de Pluvialidade

RF <- df[!is.na(df$Rainfall), ] %>% # Removing missing data
  select(Location, Rainfall, lon, lat)

aggRain <- aggregate(RF$Rainfall, by = list(Locaton=RF$Location , Longitude=RF$lon , Latitude=RF$lat), FUN=sum) # Aggregating by Location and the sum of Rainfall


YlOrBr <- c("#FFFFD4", "#FED98E", "#FE9929", "#D95F0E", "#993404")

ggmap(map) +
  stat_density_2d(
    data = aggRain,
    aes(x = Longitude, y = Latitude, fill = stat(x)),
    alpha = .1,
    bins = 30,
    geom = "polygon"
  ) +
  scale_fill_gradientn(colors = YlOrBr) +
  theme(legend.position = "left") +
  xlab("Longitude") + 
  ylab("Latitude") +
  labs(fill="Rainfall in mm")

Média de Precipitação por Mês

RD <- dataset[!is.na(dataset$Rainfall), ] %>%
  select(Date, Rainfall)

# substring(RD$Date,6,7) # Month only
RainMonth = aggregate(RD$Rainfall, by=list(Month=substring(RD$Date,6,7)), FUN=mean) # Average daily Rainfall by month
ggplot(data=RainMonth, aes(x=Month, y=x, group=1)) +
  geom_line()+
  geom_point() +
  xlab("Month") +
  ylab("Average Rainfall (mm per day)")

Criando o modelo

Removendo colunas com muitos dados faltando e depois removendo entradas incompletas.

print(colMeans(is.na(dataset))) # Percentage of missing data in each column
##          Date      Location       MinTemp       MaxTemp      Rainfall 
##    0.00000000    0.00000000    0.01020899    0.00866905    0.02241853 
##   Evaporation      Sunshine   WindGustDir WindGustSpeed    WindDir9am 
##    0.43166506    0.48009762    0.07098859    0.07055548    0.07263853 
##    WindDir3pm  WindSpeed9am  WindSpeed3pm   Humidity9am   Humidity3pm 
##    0.02906641    0.01214767    0.02105046    0.01824557    0.03098446 
##   Pressure9am   Pressure3pm      Cloud9am      Cloud3pm       Temp9am 
##    0.10356799    0.10331363    0.38421559    0.40807095    0.01214767 
##       Temp3pm     RainToday  RainTomorrow 
##    0.02481094    0.02241853    0.02245978
datasetClean <- subset(dataset, select=-c(Evaporation, Sunshine, Cloud9am, Cloud3pm, Date)) # Removing features with high NA precentage
datasetClean <- datasetClean[complete.cases(datasetClean), ] # Removing NAs

One-Hot encoding as features categóricas.

WD3<-dummy(datasetClean$WindDir3pm)
WD9<-dummy(datasetClean$WindDir9am)
WGD<-dummy(datasetClean$WindGustDir)
Loc<-dummy(datasetClean$Location)

colnames(WD3) <- paste("WindDir3pm", colnames(WD3), sep = "_")
colnames(WD9) <- paste("WindDir9am", colnames(WD9), sep = "_")
colnames(WGD) <- paste("WindGustDir", colnames(WGD), sep = "_")
colnames(Loc) <- paste("Location", colnames(Loc), sep = "_")

datasetD <- cbind(datasetClean, WD3)
datasetD <- cbind(datasetD, WD9)
datasetD <- cbind(datasetD, WGD)
datasetD <- cbind(datasetD, Loc)


datasetD$RainToday <- factor(datasetD$RainToday,
             levels = c('No', 'Yes'),
             labels = c(0, 1))

datasetD$RainTomorrow <- factor(datasetD$RainTomorrow,
                   levels = c('No', 'Yes'),
                   labels = c(0, 1))

df<-tibble(datasetD)

Separando em datasets de treino e teste.

set.seed(1337)
split = sample.split(df$RainTomorrow, SplitRatio = 0.8)
training_set = subset(df, split == TRUE)
test_set = subset(df, split == FALSE)

Treinando o classificador Random Forest.

y <- training_set$RainTomorrow
X <- subset(training_set, select=-c(RainTomorrow))

classifier = randomForest(x = X,
                          y = y,
                          ntree = 100)

classifier
## 
## Call:
##  randomForest(x = X, y = y, ntree = 100) 
##                Type of random forest: classification
##                      Number of trees: 100
## No. of variables tried at each split: 10
## 
##         OOB estimate of  error rate: 14.46%
## Confusion matrix:
##       0     1 class.error
## 0 67205  3120  0.04436545
## 1  9943 10072  0.49677742

Matriz de confusão.

y_pred = predict(classifier, newdata = subset(test_set, select=-c(RainTomorrow)))

# Confusion Matrix
cm = table(test_set$RainTomorrow, y_pred)
print(cm)
##    y_pred
##         0     1
##   0 16830   751
##   1  2413  2591

Calculando a acurácia do modelo.

# Accuracy Score
accuracy = sum(y_pred == test_set$RainTomorrow) / nrow(test_set)
print(accuracy)
## [1] 0.859907

Visualizando as curvas de ROC e Precision e Recall.

precrec_obj <- evalmod(scores =as.integer(y_pred), labels = test_set$RainTomorrow)
autoplot(precrec_obj)